import datetime as dt
import holoviews as hv
import hvplot.pandas
import hvplot.xarray
import pandas as pd
import panel as pn
import pyspedas
from pytplot import get_data
pn.extension(sizing_mode="stretch_width")
from loguru import logger
from utils import *In [1]:
THEMIS Level 2 FGM file info:
Fluxgate Magnetometer (FGM) data: - The Level 2 file includes FGE (engineering) magnetic field, FGH (high-resolution) magnetic field, FGL (low-resolution) magnetic field, and FGS (spin-resolution) magnetic field. - The FGH, FGL, and FGE data are given in Spinning Spacecraft (SSL), Despun Spacecraft (DSL), Geocentric Solar Ecliptic (GSE) and Geocentric Solar Magnetospheric (GSM) coordinates. The FGS data is given in DSL, GSE, and GSM coordinates. - Units are nanotesla.
In [38]:
def calculate_PVI_xr(vec, tau, resample_frequency=None, interval_of_averaging=None):
"""
This function calculates the Partial Variance of Increments (PVI) series for a given time series.
Parameters:
vec (np.array): The input time series, represented as a numpy array.
tau (int): The time lag, in unit `s`,typically selected to lie in the inertial range of the fluctuations.
resample_frequency (int): The resample frequency, in unit `s`. If None, defaults to tau.
interval_of_averaging (int): The number of samples over which to compute the trailing average. It's often chosen to be comparable to, or greater than, a correlation length (or time) of the signal.
Returns:
PVI_series (np.array): The resulting PVI series.
"""
# Sample the vector at the given time lag (tau)
# vec_sampled = vec.resample(time=tau).mean(dim='time')
if resample_frequency is None:
resample_frequency = tau
# Interpolate to a regular time grid
# vec_sampled = vec.resample(time=pd.Timedelta(resample_frequency, unit='s')).interpolate('linear')
vec_sampled = vec.resample(
time=pd.Timedelta(resample_frequency, unit="s")
).interpolate("linear")
# Calculate the increments in the vector
n = int(tau / resample_frequency)
increments = vec_sampled.diff(dim="time", n=n, label="lower")
# Calculate the magnitudes of these increments
mag_increments = linalg.norm(increments, dims="v_dim")
# Square the magnitudes of the increments, and compute a moving average over the specified interval
if interval_of_averaging is None:
normalized_factor = np.sqrt(np.mean(np.square(mag_increments)))
PVI_series = mag_increments / normalized_factor
else:
PVI_series = mag_increments
...
# TODO: Implement moving average
if "units" in PVI_series.attrs:
del PVI_series.attrs["units"]
PVI_series["time"] = PVI_series["time"] + pd.Timedelta(tau / 2, unit="s")
return PVI_series.rename("PVI")
def PVI_map(vec, tau_range, resample_frequency=None):
"""_summary_
Args:
vec (_type_): _description_
"""
if resample_frequency == None:
PVI_series = xr.concat(
[calculate_PVI_xr(vec, tau) for tau in tau_range], dim="tau"
)
else:
PVI_series = xr.concat(
[calculate_PVI_xr(vec, tau, resample_frequency) for tau in tau_range],
dim="tau",
)
return PVI_seriesIn [2]:
probe = "c"
start = dt.datetime(2019, 1, 7)
end = dt.date.today()
datatype = "fgs"
thx_fgs = f"th{probe}_fgs_gsm"
thx_fgl = f"th{probe}_fgl_gsm"
tstart = start.strftime("%Y-%m-%d")
tstop = (start + dt.timedelta(1)).strftime("%Y-%m-%d")
trange = [tstart, tstop]
pyspedas.themis.fgm(probe=probe, trange=trange, varnames=[thx_fgs, thx_fgl])
thx_fgs_data = get_data(thx_fgs, xarray=True)
# thx_fgl_data = get_data(thx_fgl, xarray=True)
fgl_tstart, fgl_tstop = pytplot.get_timespan(thx_fgl)
fgl_tstart = pd.Timestamp(fgl_tstart, unit="s")
fgl_tstop = pd.Timestamp(fgl_tstop, unit="s")22-Jul-23 16:17:39: Downloading remote index: http://themis.ssl.berkeley.edu/data/themis/thc/l2/fgm/2019/
22-Jul-23 16:17:39: File is current: /Users/zijin/data/themis/thc/l2/fgm/2019/thc_l2_fgm_20190107_v01.cdf
In [3]:
# Time index may contain duplicate values
thx_fgs_data = thx_fgs_data.groupby("time").first()
# Interpolate to a regular time grid
thx_fgs_data = thx_fgs_data.resample(time=pd.Timedelta(8, unit="s")).interpolate(
"linear"
)In [42]:
tau_range = [8, 16, 32, 64, 128, 256]
pvi = PVI_map(thx_fgs_data, tau_range, tau_range[0])In [43]:
pvi = pvi.assign_coords({"tau": tau_range})In [44]:
pvi.hvplot.quadmesh(x="time", y="tau", logy=True)22-Jul-23 16:41:56: /Users/zijin/mambaforge/envs/cool_solar_wind/lib/python3.10/site-packages/bokeh/core/property/bases.py:259: DeprecationWarning: elementwise comparison failed; this will raise an error in the future.
return new == old
In [2]:
date_slider = pn.widgets.DateSlider(
name="Date Slider",
start=start,
end=end,
step=1000
* 60
* 60
* 24, # (number): The selected step i the slider in milliseconds
)
overview_pane = pn.pane.PNG(
themis_image_url(start, probe, "overview_moms"), name="overview_moms"
)
orbit_pane = pn.pane.GIF(
themis_image_url(start, probe, "orbit_moon"), name="orbit_moon"
)
def update_pane_image_callback(target, event):
date = event.new
image_url = themis_image_url(date, probe, target.name)
target.object = image_url
date_slider.link(overview_pane, callbacks={"value": update_pane_image_callback})
date_slider.link(orbit_pane, callbacks={"value": update_pane_image_callback})
def index_plot(probe, date):
thx_fgm = f"th{probe}_fgs_gsm"
thx_fgl = f"th{probe}_fgl_gsm"
tstart = date.strftime("%Y-%m-%d")
tstop = (date + dt.timedelta(1)).strftime("%Y-%m-%d")
trange = [tstart, tstop]
pyspedas.themis.fgm(probe=probe, trange=trange, varnames=[thx_fgm, thx_fgl])
thx_fgm_data = get_data(thx_fgm, xarray=True)
fgl_tstart, fgl_tstop = pytplot.get_timespan(thx_fgl)
fgl_tstart = pd.Timestamp(fgl_tstart, unit="s")
fgl_tstop = pd.Timestamp(fgl_tstop, unit="s")
pvi = calculate_PVI_xr(thx_fgm_data, "10s")
mag_change = 10 * calculate_mag_change_xr(thx_fgm_data, "10s")
pvi_plot = pvi.hvplot(x="time")
threshold_line = hv.Curve([(fgl_tstart, 5), (fgl_tstop, 5)])
thx_fgm_plot = thx_fgm_data.hvplot(x="time", by="v_dim")
# thx_fgl_data = get_data(thx_fgl, xarray=True)
# thx_fgl_plot = thx_fgl_data.hvplot(x="time", by="v_dim", kind="scatter")
pvi_plot.opts(ylim=(0, 10), alpha=0.3)
threshold_line.opts(alpha=0.3)
index_plot = pvi_plot * mag_change.hvplot(x="time", alpha=0.3) * threshold_line
markdown = pn.pane.Markdown()
def update_widget(event):
temp_tstart = pd.Timestamp(event.new[0])
temp_tstop = pd.Timestamp(event.new[1])
temp_datatype = "fgs"
if temp_tstart > fgl_tstart and temp_tstop < fgl_tstop:
temp_datatype = "fgl"
markdown.object = f"""
```
{{"probe": "{probe}", "tstart": "{temp_tstart.strftime('%Y-%m-%dT%H:%M:%S')}", "tstop": "{temp_tstop.strftime('%Y-%m-%dT%H:%M:%S')}", "datatype": "{temp_datatype}"}},
```
"""
rng = hv.streams.RangeX(source=index_plot)
rng.param.watch(update_widget, "x_range")
return pn.Column(markdown, pn.Row(index_plot, thx_fgm_plot))
index_pane = hvplot.bind(lambda date: index_plot(probe, date), date_slider)
app = pn.Column(pn.Row(date_slider), index_pane, pn.Row(overview_pane, orbit_pane))
appIn [34]:
date = dt.datetime(2019, 1, 7)
tstart = date.strftime("%Y-%m-%d")
tstop = (date + dt.timedelta(1)).strftime("%Y-%m-%d")
trange = [tstart, tstop]
pyspedas.themis.fgm(probe=probe, trange=trange)29-Jun-23 12:28:55: Downloading remote index: http://themis.ssl.berkeley.edu/data/themis/thc/l2/fgm/2019/
29-Jun-23 12:28:55: File is current: /Users/zijin/data/themis/thc/l2/fgm/2019/thc_l2_fgm_20190107_v01.cdf
['thc_fgs_btotal',
'thc_fgs_gse',
'thc_fgs_gsm',
'thc_fgs_dsl',
'thc_fgl_btotal',
'thc_fgl_gse',
'thc_fgl_gsm',
'thc_fgl_dsl',
'thc_fgl_ssl',
'thc_fgh_btotal',
'thc_fgh_gse',
'thc_fgh_gsm',
'thc_fgh_dsl',
'thc_fgh_ssl']
In [3]:
import holoviews as hv
# Widgets
import pandas as pd
import panel as pn
# (index_plot + thx_fgm_plot).cols(1)
xr = thx_fgm_data
start = xr.time.data.min()
end = xr.time.data.max()
start = pd.Timestamp(start).to_pydatetime()
end = pd.Timestamp(end).to_pydatetime()
# Create a DateRangeSlider widget
date_range_slider = pn.widgets.DateRangeSlider(
name="Date Range Slider",
start=start,
end=end,
value=(start, end),
format="%H:%M:%S",
)
# Create a DatetimeRangeInput widget
datetime_range_input = pn.widgets.DatetimeRangeInput(
name="Datetime Range Input",
value=(start, end),
start=start,
end=end,
)
# Create a Markdown widget to copy
markdown = pn.pane.Markdown()
def markdown_callback(target, event):
target.object = f"""
```
event = {{
"probe": "{probe}",
"tstart": "{event.new[0]}",
"tstop": "{event.new[1]}",
"datatype": "{datatype}",
}}
```"""
# Link the values of the DateRangeSlider and DatetimeRangeInput widgets
date_range_slider.link(datetime_range_input, value="value")
date_range_slider.link(markdown, callbacks={"value": markdown_callback})
datetime_range_input.link(date_range_slider, value="value")
def plot(dt_range):
tstart = str(dt_range[0])
tstop = str(dt_range[1])
return plot_event(probe, tstart, tstop, datatype, output="all")
interactive_plot = pn.bind(plot, date_range_slider)
pn.Column(
pn.Row(date_range_slider, datetime_range_input),
markdown,
pn.Row(index_plot, thx_fgm_plot.apply.opts(xlim=date_range_slider)),
interactive_plot,
)